function [ v , zeta , Lh ] = SolveODE( Grids , Params , GE )

% Define useful functions
run MiscFunctions.m
run Define_fzeta.m

% Define empty arrays
v    = zeros(Grids.Nx,1);
zeta = zeros(Grids.Nx,1);
Lh = zeros(Grids.Nx,1);

% To make notation easier to read
DLearn = Params.DLearn ;
C  = ( 1 - Params.alpha ) * Params.psi / Params.P1 ;
E = ( 1 - Params.alpha ) * Params.P2 ;

% Initialize
v(Grids.Nx)    = GE.vU ;
zeta(Grids.Nx) = Params.zetaL ;
Lh(Grids.Nx) = Lhat( v(Grids.Nx) , zeta(Grids.Nx) , Grids.x(Grids.Nx) ) ;

% Loop
i=Grids.Nx ;
while i > 1
    
    u_i = min( max( uRate(v(i),zeta(i)) , Grids.minu ) , Grids.maxu ) ;
    G_i = G(v(i),zeta(i)) ;
    Phi_i = Phi(v(i),zeta(i)) ;
    learn0_i = Params.phi * u_i / ( DLearn + Params.phi * u_i ) / ( DLearn + Params.delta * zeta(i) ) ;
    learnA_i = u_i * Phi_i / ( Params.b + v(i) ) ;
    learnB_i = Params.delta * ( 1 - u_i ) + Phi_i * SbarPrimeOverSbar(zeta(i)) ;
    
    dzeta_i = -  Lh(i) * Theta(v(i),zeta(i)) * u_i ...
            / ( GE.C * fzeta( zeta(i) ) / Grids.fx(i) * Vacancies(v(i),zeta(i),Grids.x(i),Grids.tax(i)) )  ;
        
    dzeta_i = real(dzeta_i) ;

    zeta(i-1) = min( zeta(i) - Grids.dx(i) * dzeta_i , 2 * Params.zetaU ) ;

    % v derivative
    Ai = 1 / ( ... basic LHS term
                 Params.alpha ...
               + v(i) / ( Params.b + v(i) ) * ( zeta(i) - 1 - C ) ...
               ... migration contribution
               + v(i) * Params.eps * C * GV(v(i),zeta(i)) / G_i ...
               ... learning term
               - E * learn0_i * learnA_i ...
             ) ...
        ... basic RHS term           
        * (   ( 1 - Params.alpha ) / Grids.x(i) ...
        ... migration contribution
            + Params.eps * C * GZ(v(i),zeta(i)) / G_i * dzeta_i ...
        ... learning term
              - E * learnB_i * dzeta_i ...
        ... inefficiency term
            + Params.ineff * Params.alpha ...
              * ( SbarPrimeOverSbar( zeta(i) ) + log( Params.B0 / ( Params.b + v(i) ) ) ) ...
              * dzeta_i ...
        ... policy term
            +   ( 1 - Params.alpha ) * Grids.DlogPolicy(i) ...
          ) ;
      
    v(i-1) = min( max(  v(i) / ( 1 + Grids.vsmooth * Ai * Grids.dx(i) ) * ( 1 - ( 1 - Grids.vsmooth ) * Ai * Grids.dx(i) ) , Grids.vmin ) , Grids.vmax ) ;
    
    Lh(i-1) = Lhat( v(i-1) , zeta(i-1) , Grids.x(i-1) ) ;
                       
    i = i - 1 ;

end

end

